#include <iostream>
#include <cmath>
#include <string>
#include <fstream>

using namespace std;

const double flux = 5.22227e-10;
const double POT = 6.46165e20;
const double R = 610.6;
const double rho = 0.845;
const double deltaT = 18;
const double Na = 6.02214e23;
const double Nn = Na * rho * 4 * M_PI * R * R * R / 3;

const double ibg[51] = {0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 1.571555e-01, 1.364424e+00, 4.705236e+00, 1.758978e+01, 3.496415e+01, 6.816501e+01, 1.368511e+02, 2.066710e+02, 3.374507e+02, 5.005924e+02, 6.672170e+02, 8.545720e+02, 1.014480e+03, 1.165838e+03, 1.241060e+03, 1.300197e+03, 1.349314e+03, 1.328437e+03, 1.331230e+03, 1.305422e+03, 1.274116e+03, 1.217948e+03, 1.147513e+03, 1.085130e+03, 1.029805e+03, 1.012486e+03, 9.418484e+02, 8.868173e+02, 8.388830e+02, 7.663445e+02, 7.245175e+02, 6.952694e+02, 6.376685e+02, 6.139051e+02, 5.547155e+02, 5.363876e+02, 4.934474e+02, 4.632671e+02, 4.293046e+02, 4.191864e+02, 3.713735e+02, 3.634909e+02, 3.295256e+02, 3.028509e+02, 4.277227e+03};
const double bg[51] = {1.468132e+01,  1.751390e+02,  5.432570e+02,  6.452100e+02,  6.288804e+02,  5.888983e+02,  5.721696e+02,  5.599448e+02,  5.320706e+02,  5.089781e+02,  4.717286e+02,  4.501799e+02,  4.296122e+02,  3.620738e+02,  3.576632e+02,  3.276326e+02,  3.137620e+02,  2.857147e+02,  2.655224e+02,  2.647242e+02,  2.455104e+02,  2.285900e+02,  2.313580e+02,  2.262919e+02,  2.283181e+02,  2.124414e+02,  1.952493e+02,  1.962794e+02,  1.876680e+02,  2.057413e+02,  2.050090e+02,  1.850047e+02,  1.591487e+02,  1.803423e+02,  1.927868e+02,  2.144821e+02,  2.260172e+02,  2.518160e+02,  2.556633e+02,  2.555962e+02,  2.480172e+02,  2.373849e+02,  1.840215e+02,  1.800383e+02,  1.857559e+02,  2.038884e+02,  2.060577e+02,  2.033785e+02,  2.086383e+02,  2.217605e+02,  2.130483e+02};
const double data[51] = {4.000000e+01, 4.950000e+02, 1.902000e+03, 2.722000e+03, 2.908000e+03, 2.878000e+03, 2.934000e+03, 3.003000e+03, 2.892000e+03, 2.720000e+03, 2.651000e+03, 2.656000e+03, 2.633000e+03, 2.409000e+03, 2.364000e+03, 2.250000e+03, 2.104000e+03, 2.070000e+03, 1.929000e+03, 1.832000e+03, 1.762000e+03, 1.618000e+03, 1.546000e+03, 1.458000e+03, 1.371000e+03, 1.203000e+03, 1.174000e+03, 1.055000e+03, 9.930000e+02, 9.080000e+02, 8.320000e+02, 7.570000e+02, 6.840000e+02, 6.330000e+02, 6.520000e+02, 7.090000e+02, 7.010000e+02, 6.710000e+02, 5.890000e+02, 5.110000e+02, 4.890000e+02, 3.920000e+02, 3.590000e+02, 3.320000e+02, 3.450000e+02, 3.430000e+02, 3.560000e+02, 3.410000e+02, 3.150000e+02, 3.240000e+02, 3.360000e+02};

const double bgHEp[30] = {6.544238e+01, 6.131658e+01, 6.802268e+01, 6.614949e+01, 5.759906e+01, 5.948360e+01, 6.494799e+01, 7.288111e+01, 7.266334e+01, 8.067760e+01, 7.378685e+01, 7.559236e+01, 6.692164e+01, 5.954329e+01, 5.842109e+01, 6.112299e+01, 5.521307e+01, 5.631305e+01, 5.625658e+01, 5.765196e+01, 6.470708e+01, 6.157652e+01, 6.896612e+01, 5.997896e+01, 7.765791e+01, 6.953401e+01, 7.666214e+01, 8.304232e+01, 9.066559e+01, 1.028902e+02};
const double dataHEp[30] = {4.600000e+02, 4.570000e+02, 4.390000e+02, 4.310000e+02, 3.800000e+02, 3.830000e+02, 3.530000e+02, 3.630000e+02, 3.710000e+02, 3.960000e+02, 2.750000e+02, 2.760000e+02, 2.220000e+02, 2.030000e+02, 2.170000e+02, 2.240000e+02, 2.000000e+02, 1.840000e+02, 1.970000e+02, 1.590000e+02, 1.550000e+02, 1.620000e+02, 1.370000e+02, 1.520000e+02, 1.260000e+02, 1.180000e+02, 1.330000e+02, 1.170000e+02, 1.380000e+02, 1.230000e+02};

const double bgHEN[30] = {4.089986e+02, 3.728637e+02, 3.529947e+02, 3.474604e+02, 3.084260e+02, 2.939212e+02, 3.071228e+02, 3.304660e+02, 3.354977e+02, 3.431949e+02, 3.271199e+02, 3.149625e+02, 2.886637e+02, 2.234001e+02, 2.394022e+02, 2.461791e+02, 2.466851e+02, 2.378370e+02, 2.380596e+02, 2.461632e+02, 2.536219e+02, 2.657772e+02, 2.737242e+02, 2.553884e+02, 2.966817e+02, 2.959479e+02, 3.014583e+02, 3.316176e+02, 3.625725e+02, 3.955619e+02};
const double dataHEN[30] = {2.280000e+03, 1.972000e+03, 1.815000e+03, 1.516000e+03, 1.356000e+03, 1.163000e+03, 1.154000e+03, 1.136000e+03, 1.150000e+03, 9.980000e+02, 7.680000e+02, 7.140000e+02, 5.730000e+02, 5.010000e+02, 4.520000e+02, 4.960000e+02, 4.620000e+02, 4.150000e+02, 4.040000e+02, 4.060000e+02, 3.950000e+02, 3.930000e+02, 3.910000e+02, 4.100000e+02, 3.750000e+02, 3.760000e+02, 4.230000e+02, 4.230000e+02, 4.390000e+02, 4.380000e+02};

void read_distr(string filename, double tab[5][51]); //read true energy distribution from file
void read_distrHE(string filename, double tab[5][30]); //for HE sample
void read_rtab(string filename, double tab[51][51]); //read migration matric from file
void read_rtabHE(string filename, double tab[30][30]); //for HE sample
void read_Erec(string filename, double *tab); //read reconstructed energy distribution from file
void makeR(double *mu, double M[51][51], double R[51][51]); //create response matrix
void makeRHE(double *mu, double M[30][30], double R[30][30]); //for HE sample
void trans (double *in, double T[][51], double *out); //make reconstructed distribution from true 
void transHE (double *in, double T[][30], double *out); //for HE sample
void true2rec(double Et[5][51], double *rec); //run above
void true2recp(double Et[5][30], double *rec); //for HE sample single proton
void true2recN(double Et[5][30], double *rec); //for HE sample multiproton

